library(foreign)

library(lattice)

library(Hmisc)

setwd("/Users/rene/Dropbox/EPideal/lsq/revisions_lsq_122809")

load("posteriorsamples072009.RData")

metadata <- vector("list", 1000)
s <- vector("list", 1000)
e <- vector("list", 1000)
l <- vector("list", 1000)
v <- vector("list", 1000)

uncdata<-read.csv("UNCids.csv",header=TRUE,sep=",")
uncdata<-uncdata[1:98,1:14]

uncpos<-uncdata[,c(6,11)]

for (i in 1:1000) {

metadata[[i]] <- merge(z[[i]],uncpos,by.x="NP")

}

for (i in 1:1000) {
 yearunc<-metadata[[i]][,11]*metadata[[i]][,2]
 year2unc<-metadata[[i]][,11]*metadata[[i]][,9]
 logodds<-1-log((10/metadata[[i]][,10])-1)
 metadata[[i]]<-cbind(metadata[[i]],yearunc,year2unc,logodds)
}

for (i in 1:1000) {

s[[i]]<-metadata[[i]][metadata[[i]][,6]=="S",]
e[[i]]<-metadata[[i]][metadata[[i]][,6]=="E",]
l[[i]]<-metadata[[i]][metadata[[i]][,6]=="L",]
v[[i]]<-metadata[[i]][metadata[[i]][,6]=="V",]

}


################################################################################


#summary(s[[1]]$position)


year <- seq(1, 5, 0.001)

position <- rep(4.920, 4001)

s.new.low <- as.data.frame(cbind(year,position))

s.new.low$year2 <- (s.new.low$year)^2

s.new.low$yearunc <- s.new.low$year*s.new.low$position

s.new.low$year2unc <- s.new.low$year2*s.new.low$position


year <- seq(1, 5, 0.001)

position <- rep(6.860, 4001)

s.new.high <- as.data.frame(cbind(year,position))

s.new.high$year2 <- (s.new.high$year)^2

s.new.high$yearunc <- s.new.high$year*s.new.high$position

s.new.high$year2unc <- s.new.high$year2*s.new.high$position


################################################################################

s.r <- s[runif(100,1,1000)]

pred.s.low <- vector("list", 100)

pred.s.high <- vector("list", 100)

s.lm <- vector("list", 100)

s.short <- vector("list", 100)

for (i in 1:100) {

s.short[[i]] <- s.r[[i]][,c(10,14,2,11,9,12,13)]

names(s.short[[i]]) <- c("estimate", "logodds", "year", "position", "year2", "yearunc", "year2unc")

s.lm[[i]] <- lm(estimate ~ year + position + year2 + yearunc + year2unc, data=s.short[[i]])

pred.s.low[[i]] <- predict(s.lm[[i]], s.new.low, interval="confidence")

pred.s.high[[i]] <- predict(s.lm[[i]], s.new.high, interval="confidence")

}

ind <- vector("list", 100)

year <- vector("list", 100)

for (i in 1:100) {

ind[[i]] <- rep(i,4001)

year[[i]] <- seq(1, 5, 0.001)

}


coef.credible.s <- matrix(NA, 100, 6)

for (i in 1:100) {
 coef.credible.s[i,] <- s.lm[[i]]$coefficients
}

colnames(coef.credible.s) <- names(s.lm[[1]]$coefficients)

credible.s <- t(apply(coef.credible.s, 2, quantile, c(0.025, 0.5, 0.975)))


random.select <- round(runif(1,1,100), digits=0)

#pdf("residuals_s.pdf")

#par(mfrow=c(1,2))

#plot(s.lm[[random.select]]$fit,s.lm[[random.select]]$res,xlab="Fitted",ylab="Residuals")

#qqnorm(s.lm[[random.select]]$res,ylab="Raw Residuals")

#qqline(s.lm[[random.select]]$res)

#dev.off()


s.low <- do.call(rbind, pred.s.low)

s.low <- as.data.frame(s.low)

s.low$ind <- do.call(c, ind)

s.low$year <- do.call(c, year)

s.low$group <- rep(1, nrow(s.low))


s.high <- do.call(rbind, pred.s.high)

s.high <- as.data.frame(s.high)

s.high$ind <- do.call(c, ind)

s.high$year <- do.call(c, year)

s.high$group <- rep(2, nrow(s.high))


s.all <- rbind(s.low, s.high)

s.all$fit.t <- 10*exp(s.all$fit)/(exp(1)+exp(s.all$fit))

s.all$lwr.t <- 10*exp(s.all$lwr)/(exp(1)+exp(s.all$lwr))

s.all$upr.t <- 10*exp(s.all$upr)/(exp(1)+exp(s.all$upr))


socialists <- s.all[s.all$ind==random.select,]

socialists$pg <- rep(3, nrow(socialists))


#pdf("socialists.pdf")

xYplot(Cbind(fit.t,lwr.t,upr.t) ~ year | as.factor(ind), group=as.factor(group), type="l", method="bands", label.curves=F, col=c("black", "grey51"), ylab="Fit", xlab="Year", par.settings=list(strip.background=list(col=c(0))), strip=strip.custom(strip.names=TRUE, var.name="Draw"), data=s.all[s.all$ind==1,])

#dev.off()


################################################################################


#summary(e[[1]]$position)


year <- seq(1, 5, 0.001)

position <- rep(2, 4001)

s.new.low <- as.data.frame(cbind(year,position))

s.new.low$year2 <- (s.new.low$year)^2

s.new.low$yearunc <- s.new.low$year*s.new.low$position

s.new.low$year2unc <- s.new.low$year2*s.new.low$position


year <- seq(1, 5, 0.001)

position <- rep(6.88, 4001)

s.new.high <- as.data.frame(cbind(year,position))

s.new.high$year2 <- (s.new.high$year)^2

s.new.high$yearunc <- s.new.high$year*s.new.high$position

s.new.high$year2unc <- s.new.high$year2*s.new.high$position


#######


s.r <- e[runif(100,1,1000)]

pred.s.low <- vector("list", 100)

pred.s.high <- vector("list", 100)

s.lm <- vector("list", 100)

s.short <- vector("list", 100)

for (i in 1:100) {

s.short[[i]] <- s.r[[i]][,c(10,2,11,9,12,13)]

names(s.short[[i]]) <- c("estimate", "year", "position", "year2", "yearunc", "year2unc")

s.lm[[i]] <- lm(estimate ~ year + position + year2 + yearunc + year2unc, data=s.short[[i]])

pred.s.low[[i]] <- predict(s.lm[[i]], s.new.low, interval="confidence")

pred.s.high[[i]] <- predict(s.lm[[i]], s.new.high, interval="confidence")

}


ind <- vector("list", 100)

year <- vector("list", 100)

for (i in 1:100) {

ind[[i]] <- rep(i,4001)

year[[i]] <- seq(1, 5, 0.001)

}


coef.credible.e <- matrix(NA, 100, 6)

for (i in 1:100) {
 coef.credible.e[i,] <- s.lm[[i]]$coefficients
}

colnames(coef.credible.e) <- names(s.lm[[1]]$coefficients)

credible.e <- t(apply(coef.credible.e, 2, quantile, c(0.025, 0.5, 0.975)))


random.select <- round(runif(1,1,100), digits=0)

pdf("residuals_e.pdf")

par(mfrow=c(1,2))

plot(s.lm[[random.select]]$fit,s.lm[[random.select]]$res,xlab="Fitted",ylab="Residuals")

qqnorm(s.lm[[random.select]]$res,ylab="Raw Residuals")

qqline(s.lm[[random.select]]$res)

dev.off()


s.low <- do.call(rbind, pred.s.low)

s.low <- as.data.frame(s.low)

s.low$ind <- do.call(c, ind)

s.low$year <- do.call(c, year)

s.low$group <- rep(1, nrow(s.low))


s.high <- do.call(rbind, pred.s.high)

s.high <- as.data.frame(s.high)

s.high$ind <- do.call(c, ind)

s.high$year <- do.call(c, year)

s.high$group <- rep(2, nrow(s.high))


s.all <- rbind(s.low, s.high)


conservatives <- s.all[s.all$ind==random.select,]

conservatives$pg <- rep(4, nrow(conservatives))


pdf("conservatives.pdf")

xYplot(Cbind(fit,lwr,upr) ~ year, group=as.factor(group), type="l", method="bands", label.curves=F, col=c("black", "grey51"), ylab="Fit", xlab="Year", par.settings=list(strip.background=list(col=c(0))), strip=strip.custom(strip.names=TRUE, var.name="Draw"), data=s.all[s.all$ind==1,])

dev.off()


################################################################################


#summary(l[[1]]$position)


year <- seq(1, 5, 0.001)

position <- rep(2.31, 4001)

s.new.low <- as.data.frame(cbind(year,position))

s.new.low$year2 <- (s.new.low$year)^2

s.new.low$yearunc <- s.new.low$year*s.new.low$position

s.new.low$year2unc <- s.new.low$year2*s.new.low$position


year <- seq(1, 5, 0.001)

position <- rep(7, 4001)

s.new.high <- as.data.frame(cbind(year,position))

s.new.high$year2 <- (s.new.high$year)^2

s.new.high$yearunc <- s.new.high$year*s.new.high$position

s.new.high$year2unc <- s.new.high$year2*s.new.high$position


#######


s.r <- l[runif(100,1,1000)]

pred.s.low <- vector("list", 100)

pred.s.high <- vector("list", 100)

s.lm <- vector("list", 100)

s.short <- vector("list", 100)

for (i in 1:100) {

s.short[[i]] <- s.r[[i]][,c(10,2,11,9,12,13)]

names(s.short[[i]]) <- c("estimate", "year", "position", "year2", "yearunc", "year2unc")

s.lm[[i]] <- lm(estimate ~ year + position + year2 + yearunc + year2unc, data=s.short[[i]])

pred.s.low[[i]] <- predict(s.lm[[i]], s.new.low, interval="confidence")

pred.s.high[[i]] <- predict(s.lm[[i]], s.new.high, interval="confidence")

}


ind <- vector("list", 100)

year <- vector("list", 100)

for (i in 1:100) {

ind[[i]] <- rep(i,4001)

year[[i]] <- seq(1, 5, 0.001)

}


coef.credible.l <- matrix(NA, 100, 6)

for (i in 1:100) {
 coef.credible.l[i,] <- s.lm[[i]]$coefficients
}

colnames(coef.credible.l) <- names(s.lm[[1]]$coefficients)

credible.l <- t(apply(coef.credible.l, 2, quantile, c(0.025, 0.5, 0.975)))


random.select <- round(runif(1,1,100), digits=0)

pdf("residuals_l.pdf")

par(mfrow=c(1,2))

plot(s.lm[[random.select]]$fit,s.lm[[random.select]]$res,xlab="Fitted",ylab="Residuals")

qqnorm(s.lm[[random.select]]$res,ylab="Raw Residuals")

qqline(s.lm[[random.select]]$res)

dev.off()


s.low <- do.call(rbind, pred.s.low)

s.low <- as.data.frame(s.low)

s.low$ind <- do.call(c, ind)

s.low$year <- do.call(c, year)

s.low$group <- rep(1, nrow(s.low))


s.high <- do.call(rbind, pred.s.high)

s.high <- as.data.frame(s.high)

s.high$ind <- do.call(c, ind)

s.high$year <- do.call(c, year)

s.high$group <- rep(2, nrow(s.high))


s.all <- rbind(s.low, s.high)


liberals <- s.all[s.all$ind==random.select,]

liberals$pg <- rep(1, nrow(liberals))


pdf("liberals.pdf")

xYplot(Cbind(fit,lwr,upr) ~ year | as.factor(ind), group=as.factor(group), type="l", method="bands", label.curves=F, col=c("black", "grey51"), ylab="Fit", xlab="Year", par.settings=list(strip.background=list(col=c(0))), strip=strip.custom(strip.names=TRUE, var.name="Draw"), data=s.all[s.all$ind==1,])

dev.off()


################################################################################


#summary(v[[1]]$position)


year <- seq(1, 5, 0.001)

position <- rep(1.83, 4001)

s.new.low <- as.data.frame(cbind(year,position))

s.new.low$year2 <- (s.new.low$year)^2

s.new.low$yearunc <- s.new.low$year*s.new.low$position

s.new.low$year2unc <- s.new.low$year2*s.new.low$position


year <- seq(1, 5, 0.001)

position <- rep(6.36, 4001)

s.new.high <- as.data.frame(cbind(year,position))

s.new.high$year2 <- (s.new.high$year)^2

s.new.high$yearunc <- s.new.high$year*s.new.high$position

s.new.high$year2unc <- s.new.high$year2*s.new.high$position


#######


s.r <- v[runif(100,1,1000)]

pred.s.low <- vector("list", 100)

pred.s.high <- vector("list", 100)

s.lm <- vector("list", 100)

s.short <- vector("list", 100)

for (i in 1:100) {

s.short[[i]] <- s.r[[i]][,c(10,2,11,9,12,13)]

names(s.short[[i]]) <- c("estimate", "year", "position", "year2", "yearunc", "year2unc")

s.lm[[i]] <- lm(estimate ~ year + position + year2 + yearunc + year2unc, data=s.short[[i]])

pred.s.low[[i]] <- predict(s.lm[[i]], s.new.low, interval="confidence")

pred.s.high[[i]] <- predict(s.lm[[i]], s.new.high, interval="confidence")

}


ind <- vector("list", 100)

year <- vector("list", 100)

for (i in 1:100) {

ind[[i]] <- rep(i,4001)

year[[i]] <- seq(1, 5, 0.001)

}


coef.credible.g <- matrix(NA, 100, 6)

for (i in 1:100) {
 coef.credible.g[i,] <- s.lm[[i]]$coefficients
}

colnames(coef.credible.g) <- names(s.lm[[1]]$coefficients)

credible.g <- t(apply(coef.credible.g, 2, quantile, c(0.025, 0.5, 0.975)))


random.select <- round(runif(1,1,100), digits=0)

pdf("residuals_g.pdf")

par(mfrow=c(1,2))

plot(s.lm[[random.select]]$fit,s.lm[[random.select]]$res,xlab="Fitted",ylab="Residuals")

qqnorm(s.lm[[random.select]]$res,ylab="Raw Residuals")

qqline(s.lm[[random.select]]$res)

dev.off()


s.low <- do.call(rbind, pred.s.low)

s.low <- as.data.frame(s.low)

s.low$ind <- do.call(c, ind)

s.low$year <- do.call(c, year)

s.low$group <- rep(1, nrow(s.low))


s.high <- do.call(rbind, pred.s.high)

s.high <- as.data.frame(s.high)

s.high$ind <- do.call(c, ind)

s.high$year <- do.call(c, year)

s.high$group <- rep(2, nrow(s.high))


s.all <- rbind(s.low, s.high)


greens <- s.all[s.all$ind==random.select,]

greens$pg <- rep(2, nrow(greens))


pdf("greens.pdf")

xYplot(Cbind(fit,lwr,upr) ~ year | as.factor(ind), group=as.factor(group), type="l", method="bands", label.curves=F, col=c("black", "grey51"), ylab="Fit", xlab="Year", par.settings=list(strip.background=list(col=c(0))), strip=strip.custom(strip.names=TRUE, var.name="Draw"), data=s.all[s.all$ind==1,])

dev.off()


################################################################################


allparties <- rbind(socialists,conservatives,liberals,greens)

levels <- c("Liberals","Greens","Socialists","Conservatives")

allparties$pgf <- factor(allparties$pg, labels = levels)


pdf("predall.pdf")

xYplot(Cbind(fit,lwr,upr) ~ year | pgf, group=as.factor(group), type="l", method="bands", label.curves=F, col=c("black", "grey51"), ylab="Fit", xlab="Year", par.settings=list(strip.background=list(col=c(0))), strip=strip.custom(strip.names=TRUE, var.name="EPG"), data=allparties)

dev.off()


################################################################################


postsum <- function(x,n) {

m <- rowMeans(x[,-c(ncol(x):(ncol(x)-(ncol(x)-n-1)))])

q <- apply(x[,c(1:n)],1,quantile, probs = c(0.05, 0.95)) 

z <- cbind(x[,-c(1:n)],q[1,],m,q[2,])

colnames(z)[(ncol(z)-2):(ncol(z))] <- c("lower","estimate","upper")

z

}


sum.c <- postsum(master.c,10000)

partyleader <- idata[,c(1,9)]

sum.c <- merge(sum.c, partyleader, by.x="MEPID")

c.sds <- aggregate(sum.c$estimate,list(MEPID=sum.c$MEPID),sd)

colnames(c.sds)[2] <- "sdest"

sum.c <- merge(sum.c, c.sds, by.x="MEPID")

#summary(lm(sdest ~ partyleader, data=sum.c[sum.c$year==1,]))

t.test(sum.c[sum.c$partyleader==0 & sum.c$year==1,12],sum.c[sum.c$partyleader==1 & sum.c$year==1,12], alternative="greater")




sum.s <- postsum(master.s,10000)

sum.s <- merge(sum.s, partyleader, by.x="MEPID")

s.sds <- aggregate(sum.s$estimate,list(MEPID=sum.s$MEPID),sd)

colnames(s.sds)[2] <- "sdest"

sum.s <- merge(sum.s, s.sds, by.x="MEPID")

#summary(lm(sdest ~ partyleader, data=sum.s[sum.s$year==1,]))

t.test(sum.s[sum.s$partyleader==0 & sum.s$year==1,12],sum.s[sum.s$partyleader==1 & sum.s$year==1,12], alternative="greater")



sum.l <- postsum(master.l,10000)

sum.l <- merge(sum.l, partyleader, by.x="MEPID")

l.sds <- aggregate(sum.l$estimate,list(MEPID=sum.l$MEPID),sd)

colnames(l.sds)[2] <- "sdest"

sum.l <- merge(sum.l, l.sds, by.x="MEPID")

#summary(lm(sdest ~ partyleader, data=sum.l[sum.l$year==1,]))

t.test(sum.l[sum.l$partyleader==0 & sum.l$year==1,12],sum.l[sum.l$partyleader==1 & sum.l$year==1,12], alternative="greater")



sum.g <- postsum(master.g,10000)

sum.g <- merge(sum.g, partyleader, by.x="MEPID")

g.sds <- aggregate(sum.g$estimate,list(MEPID=sum.g$MEPID),sd)

colnames(g.sds)[2] <- "sdest"

sum.g <- merge(sum.g, g.sds, by.x="MEPID")

#summary(lm(sdest ~ partyleader, data=sum.g[sum.g$year==1,]))

t.test(sum.g[sum.g$partyleader==0 & sum.g$year==1,12],sum.g[sum.g$partyleader==1 & sum.g$year==1,12], alternative="greater")

